All nonspherical perturbations of the Choptuik spacetime decay 
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V • I. INTRODUCTION AND SUMMARY 

in 

We have many and powerful results about the static or stationary end-states of gravitational collapse. However, 
very little is known in comparison about the dynamical evolution towards them. Analytical studies are limited by the 
complicated non-linear nature of the equations. Numerical studies can fill this gap if they can demonstrate generic 
. behavior. 

^\ ' Starting with the pioneering work of Choptuik jlj], a number of authors have shown that, despite the complicated 
nature of the equations, the threshold of gravitational collapse is strikingly simple (|]. Following the initial ideas of 
Evans || it has been possible to explain this simplicity as the consequence of the existence of a "critical solution" 
which acts as an intermediate attractor in phase space. This solution has a single linearly-unstable eigenmode which 
drives out every nearby solution cither towards black hole formation or dispersal, leaving flat space behind. 

This body of work expands our understanding of the dynamical process of collapse, borrowing concepts and tools 
from the theory of dynamical systems. The emphasis is shifted to phase space, and within it, to solutions with special 
stability characteristics: 



cr 

$_i , 1. First, we look for global attractors. The Minkowski and Kerr-Newman solutions are the only possible end-states 

of collapse in the Einstein-scalar-Maxwell system. 



2. Then we look for codimension-one attractors, which separate phase space into basins of attraction of the global 
attractors. These solutions are also very important. For example, the study of the trajectories connecting 
the codimension-one attractors with the global attractors gives us a qualitative picture of marginal collapse 
because many different trajectories tend to approach them and arrive at the attractors along them. In this 
terminology, Choptuik discovered the first codimension-one attractor. For reasons still unknown, many of the 
codimension-one attractors are self-similar. 

3. The long-term objective is the construction of a picture of the unfolding of trajectories in phase space. It would 
contain all the dynamical information about a given system. Furthermore it is the natural place to accommodate 
the zoo of special solutions we currently know of, including naked singularities. 

In a previous paper [|| we addressed the question of whether the Choptuik solution was a codimension-one solution 
in the system Einstein-Maxwell-charged scalar field, restricted to spherical symmetry, and obtained an affirmative 
answer, which has been confirmed in independent work (5). In this paper we address the same question for the system 
Einstein-real scalar field, but this time allowing for arbitrary small deviations from spherical symmetry, and we 
obtain again an affirmative answer. The study of the Einstein-Maxwell-charged scalar field system beyond spherical 
symmetry will be reported elsewhere. 

This result, together with a parallel result on the collapse of a perfect fluid ||, strongly suggests that critical 
phenomena in gravitational collapse are still present in the absence of spherical symmetry. An equally strong indication 
that critical phenomena are not restricted to spherical symmetry is provided by numerical work on the critical collapse 
of axisymmetric vacuum spacetimes , which shows universality and scaling similar to that of the spherical scalar 
field. 
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The plan of the article is as follows: In section II we give a complete review of the Gerlach and Sengupta 




formalism of gauge-invariant perturbations around a general spherically symmetric spacetime (which typically contains 
matter and is time-dependent). In section III we re-express these still general tensor equations of section II in an 
arbitrary basis to facilitate the study of their causal structure. The equations in this section will be of help in any study 
of linear perturbations around spherical symmetry for arbitrary matter content, in an arbitrary background coordinate 
system. In section IV we specialize the formalism to massless scalar field matter. The background solution is briefly 
reviewed in section V, where we specialize to a particular basis, and choose a coordinate system. In sections VI and 
VII we split the odd and even linear perturbation equations, respectively, into evolution equations and constraints, 
and identify free data. Section VIII describes our numerical results in detail. The appendix contains a description of 
our numerical methods for computing the background and then the perturbations on it. 

To summarize our main result here, all non-spherical physical perturbations of Choptuik's solution decay, and 
therefore the critical phenomena at the black hole threshold in scalar field collapse — universality, echoing and scaling 
- are expected to persist for initial data that deviate (slightly) from spherical symmetry. Nevertheless, the / = 2 
even-parity perturbations decay quite slowly, and may become apparent in non-spherical collapse situations. 



II. REVIEW OF GERLACH AND SENGUPTA FORMALISM OF GAUGE-INVARIANT 

PERTURBATIONS 

In this section we give a brief introduction to the formalism of Gerlach and Sengupta for perturbations around 
the most general spherically symmetric spacetime. Spacetime is decomposed as M 4 = M 2 x S 2 , where S 2 is the 
two-sphere and M 2 is a two-dimensional manifold with boundary. Tensor indices on M 4 are Greek letters, tensor 
indices on M 2 are upper case Latin letters, and tensor indices on S 2 are lower case Latin letters. We write the general 
spherically symmetric metric as 

g MV da?dx" = g AB {x D )dx A dx B + r 2 (x D )j ab (x d )dx a dx b , (1) 

where gAB is a metric and r is a scalar field on M 2 . ^ a b is the unit Gaussian-curvature metric on 5 2 . r = identifies 
the center of the spherical symmetry, where each S* 2 degenerates to a point, r — is the boundary of M 2 . In the 
same way we decompose the spherically symmetric stress-energy tensor: 

t^dx^dx" = t AB (x D )dx A dx B + Q(x D )r 2 (x D )jab(x d )dx a dx b . (2) 

For simplifying the field equations, it is useful to introduce a vector and a scalar on M 2 derived from the scalar r: 

va = , (3) 

r 

V = -r~ 2 + 2v A lA + 3v A v A . (4) 

We distinguish covariant derivatives on M 4 , M 2 and S 2 : 

g^v-x = 0, g A B\c = 0, lab:c = 0. (5) 

We shall also need the covariantly constant unit antisymmetric tensors with respect to gAB and g a b, which we call 
e A B and e ab . 

The Einstein equations in spherical symmetry are 

Gab = ~2(v A \b + v A v B ) + gAsVa = 8irt A B, (6) 
~G a a = -n + VA lA +v A v A =8ttQ, (7) 

where G a a denotes the partial trace over G p „. TZ = \R^a A IS ^ ne Gaussian-curvature scalar of M 2 . The four- 
dimensional Ricci scalar is R = 2(1Z — Vo). The conservation equation for the stress-energy tensor in spherical 
symmetry is 

t A B lB + 2t AB v B =2Qv A . (8) 

As a manifestation of the contracted Bianchi identities, (Q) can be obtained as a derivative of equation (|^) , provided 
that (§) holds. 

Now we introduce an arbitrary (not spherically symmetric) perturbation of this spacetime: — * g^v{x D ) + 
hfj,^(x D , x d ) and again we perform a 2+2 decomposition. Furthermore we decompose the angular (x d ) dependence 
into series of tensorial spherical harmonics: 
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• Y l m (x d ) are the scalar spherical harmonics, 

• the objects Y™, and Sf 1 = e a h ^ m . b form a complete basis of vector harmonics, and 

• following Zerilli JO], we use the following basis of symmetric tensor harmonics: Yj m 7 i,, Z™ afc = ^ m . Q (, + 
K'+^ ym^^ an j S™ a . b + S'[ n b . al which is a linear combination of the basis introduced by Regge and Wheeler 
p2|. For I = 0, 1 there is only one linearly independent tensor, namely 7 a f,YJ m , while the other two tensors 
vanish. Gerlach and Sengupta initially Ji| used the Regge- Wheeler basis, but in || changed to Zerilli's basis in 
order to include the cases I = 0, 1 into a single formalism. 

All these spherical harmonics have definite parity under spatial inversion: a spherical harmonic with label I is called 
even if it has parity (—1)' and odd if its parity is (— Y™, Y, m ' and ZY 1 ■ are even and ST 1 and S™ b + SJ n b . a 

are odd. (An alternative terminology is polar instead of even, and axial instead of odd.) Even and odd perturbations 
decouple, and different values of I and to decouple. Furthermore, the perturbation equations do not depend on m. 
In the following we consider one value of I and to at a time, and suppress both the indices I and m and the explicit 
summation over them, is decomposed into 

tiAB = h AB Y, (9) 
h Ab = h%Y b + h%S b , (10) 
h ab = r 2 K lab Y + r 2 GZ ab + h{S a:b + S bta ). (11) 

Note that the left-hand sides are components of a tensor on M 4 . On the right-hand side Hab is a tensor on M 2 , and 
Y is a scalar on S 2 . Similar remarks apply to the other definitions. In the same way we decompose the perturbation 
At^y into tensorial spherical harmonics: 

At AB = At AB Y, (12) 
At Ab = At^Y b + At^S b , (13) 
At ab = r 2 At 3 lab Y + At 2 Z ab + At(S a:b + S b:a ). (14) 

where we use the superindices 2 and 3 in order to follow the notation of J^]. (They are just labels, not components of 
any vector.) Some of the coefficients on the right hand side of these expansions are not defined for / = 0, 1 because 
the corresponding spherical harmonics vanish. In the following, we always point out which of the general equations 
continue to hold for / = and / = 1 if one sets these coefficients to zero. 

Now we define gauge- invariant variables, which do not contain perturbations of the background generated by simple 
coordinate transformations on this background. With the shorthand 

PA = h^-^r 2 G {A , (15) 
a complete set of gauge-invariant metric perturbations is 

kAB = h A B - (jpA\B +Pb\a), (16) 

k A = h ( i-h\ A + 2hv A , (17) 

k = K+ l l±}l G -2v A pA. (18) 
A complete set of gauge-invariant matter perturbations is 

~ tAB\CP C - tACP C \ B - tBCP C \ A , (19) 

t AC p C (20) 

^(r 2 Q) l cP C + l -^ 1 QG, (21) 

r 2 QG, (22) 
Qh A , (23) 
}h. (24) 



Tab = 


Ai A B 


T A = 


Atl- 


rp3 


At 3 - 


rpl 


At 2 ~ 


L A = 


At°- 


L = 


At - i 



3 



(These are only partially gauge-invariant for I — 0, 1, and therefore a partial gauge fixing is required in those cases in 
order to eliminate coordinate transformations from the set of arbitrary perturbations.) 

For I > 2 the metric and matter perturbations in any particular gauge can be obtained by freely choosing values 
for G, and h (the gauge in which these all vanish is Regge- Wheeler gauge jl^]) and then solving the definitions of 
the gauge-invariant perturbations for the "naked" perturbations. 

The perturbed Einstein equations, expressed only in gauge-invariant perturbations, are 



I > 



(kcA\B + k C B\A - k AB \ C )v c - g AB (2k CD lD - k D D i c )v 



c 



n , ^ fir K l + ^)\ 

-{k\ A v B + k\ B v A + k\ AB ) + I Vq H I 



k AB 



9AB 



k F ^ 2 ~ 



2k DF v DlF + 3k DF v D v F - k\ 



\{-k AU \ AB + k A A lB lB - 2k AM \ A v B + k A A]B v B + R AB (k AB - kg AB ) 



3h F v f + 



2r 2 



8irT A1 



2r 2 A ' K \A 



2k\ A v A } = 8ttT 3 



(25) 



(26) 



I > 1 



1 

2^2 



^{ k AB^ 



- k B B \ A + k B B v A - h A ) = 8irT A , 



\c 



(I- l)(/ + 2) 
2r 2 



k A = 8nL A , 



(27) 
(28) 



I > 2 



5*V 



8ttT 2 , 

8ttZ. 



(29) 
(30) 



R in ([26|) are the AB components of the four-dimensional Ricci tensor. Equations (|26j) and ( |30| ) can be obtained 
as derivatives of the other equations using the linearized equations of stress-energy conservation, which are 

I > : 



y^ 2 T Al 



1(1 + 1) 



T A - 2v A T 3 = (t AB \ D + 2t AB v D )k BD + Q(k\ A - 2kv A ) - t AB k^ 



1 BCl 1 
--jt k BC \ A - —t AB kp 



F\B 



tABk BF \ F: 



(31) 



I > 1 



r z 2 r z 2 

(r 2 L B )\ B = (l-l)(l + 2)L. 



Qik- 



-k A A 



(32) 
(33) 



III. PERTURBATION EQUATIONS FOR ARBITRARY MATTER IN AN ARBITRARY 

ORTHONORMAL BASIS 

Both in order to transform tensor equations into sets of scalar equations, and in order to separate evolution equations 
from constraints, it is desirable to introduce an orthonormal frame in M 2 , namely 

- u A u A = n A n A ee 1, u A n A = 0. (34) 
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In the presence of curvature, this cannot be a coordinate basis: 

[n, u] A = (in A — vu A , (i = u A \ A , v = n A \ A . (35) 

We define an associated basis of 2-tensors, 

9AB = -u A u B +n A n Bl e AB = n A u B - u A n Bl (36) 
p AB = u A u B + n A n B , Qab = n A u B + u A n B . (37) 

and use it to decompose the gauge-invariant metric perturbation: 

kAB = V9A B + 4>PA B + ipqAB- (38) 

We define derivatives along the basis vectors: 

/ = u A f ]Al f = n A f lAl (39) 

and re-express the even-perturbation equations in this basis. We also introduce the notation 

u A v A = U, n A v A = W, W 2 -U 2 = v A v A = v 2 . (40) 

For reference we give the background Einstein equations in frame components: 

W - U + vW - (iU - 2U 2 + 2W 2 - r- 2 = 8ir^t A A , (41) 

-W -U + vW + (iU -U 2 -W 2 = Sn^ PAB t AB , (42) 

-U' - W + (iW + vU - 2UW = &TT^q AB t AB , (43) 

-K + W - U + vW - (iU - U 2 + W 2 = 8itQ . (44) 

We use them among other things to bring all perturbation equations into a standard form by eliminating the derivatives 
of U and W. 

The complete Einstein equations for the even perturbations, still for arbitrary matter content, expressed in gauge- 
invariant variables, and decomposed in an arbitrary frame, are 

I > : 

-k + k" + vk' -nk- 2U[2k + + $ + 2(i<P + 2vip] - 2W[-2k' + <f>' + ip + 2v<p + 2(1^] 

- ^-^- (k + V) + ^(k -V)~ MU 2 + W 2 ) - A^UW + 167r(#AB + HA B )t AB = &irT A A , (45) 
-k - k" + vk' + (ik + 2U[-k + f) + tp' + 2(jL<p] + 2W[-k' + r)'-tp- 2v<p\ 

+20 (v + l ^r-) = 8np AB T AB , (46) 
— (k)' - (k'J + nk' + vk + 2U[-k' + r/ -</>'- 2(ii\}\ + 2W[-k + f) + + 2viP] 

-^( V o+ l -^fi) =8nq AB T AB , (47) 
-k - rj - + k" + rj' - 0" - (ip)' - (ip 1 ) + vk! - (j,k + vr{ - (if) - 3(v<fi' + vip + (lip' + (i(p) 
-2t/(0 + k + tp') — 2W(<p' -k' + ip)- 2(v'4> + vip + (i'tp + (up) - l \ 

-2(k - n)[n - W + u - vW + fiu + u 2 - w 2 } 

-2<P[W' + U + vW + (iU + U 2 + W 2 + (i 2 + v 2 } 

-2ip[U' + W + (jW + vU + 2UW + 2(iv] = 167rT 3 , (48) 

l>l: k + i) + <j) + iP' + 2(i(f) + 2viP - 2Ur] = -16nu A T A , (49) 

k' + T)' - (p' - ip - 2v<P - 2(iiP - 2Wf) = -16im A T Al (50) 

I > 2 : r)= -8ttT 2 . (51) 
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We have now turned the even perturbation equations into scalar form. They are already in first-order form if one 
counts /, / and /' as separate variables linked by certain trivial equations. In the following we always imply this 
first-order interpretation. The final step of the analysis is to separate the equations into evolution equations and 
constraints. This cannot be done in general, as the causal structure of the equations depends on the matter content. 

Attention must be paid to the regularity of the perturbations at r = 0. Changing to Cartesian coordinates one can 
see that regular perturbations scale as 



X 



= 6-k- 



r l rj, 
r l k, 



V 
k 

1> 



(52) 
(53) 
(54) 
(55) 



where the barred variables are 0(1) at the center. 

The odd metric-perturbations are contained in k A - We can transform the vector equation (^8|) into a scalar equation 
using the curl of k A : 



n 



JB/ -2 



(r- 2 k A )\ 



(56) 



It is possible to reconstruct k A from II for I > 2 using equation fl28|) . Therefore II alone characterizes the physical 
odd metric perturbations. For I > 2 it obeys the "odd parity master equation" M 



I > 2 



i(^n) 1 



(l-l)(Z + 2) 



n 



8ttc ab L 



MB- 



(57) 



\A 



This equation is a generalization of the Regge- Wheeler equation jl2). If we define the object Hrw = r3 n then the 
master equation is: 



n w |A A+ [V A \A 



V VA 5 



n 



-16irre AB L A{B 



(58) 



which, for Schwarzschild background in radial coordinates, and using the "tortoise" coordinate r* , is the Regge- 
Wheeler equation: 



d 2 U 



d 2 n 



dt 2 dr* 2 
We enforce regularity at the origin by defining 



2M\ (6M + 
J- — 3 5 — y^RW 

r J \ r° r J 



= 0. 



n = r /_2 n, 



(59) 



(60) 



and equation ( p7| ) in an arbitrary basis becomes 

i > 2 : -ft + n" + \v + 2(1 + l)W]n' - \p + 2(1 + l)C/]ft 



+(/ + 2) 



+ (I - l)(v A v A - r- 2 ) 



n = -16Trr-'[-{u A L A )' + (n A L A ) + [in A L A - vu A L A \. (61) 



In the special case I = 1, k A is defined by II only up to a gradient, but precisely this gradient is a gauge degree of 
freedom, so that II again contains all the gauge-invariant information. As we have (r 2 L A )\ A = 0, L A can be expressed 
as r 2 L A = e AB T\ B , with T a new scalar. Equation d2q) can be integrated to obtain the algebraic Einstein equation 



/ = 1 



r 4 n = 16ttT + const. 



(62) 



The integration constant must be zero if the perturbed spacetime is to be regular in r = 0. (If the background 
spacetime is Schwarzschild, then this integration constant parameterizes an infinitesimal angular momentum taking 
Schwarzschild into Kerr.) For I — 0, there are no odd-parity perturbations at all. 
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IV. THE MASSLESS SCALAR FIELD MODEL 



In the remainder of the paper, we restrict attention to a particular matter model, the real massless scalar field tp 
with stress-energy tensor 

tuv = V,nV,u - \aiiv9,x<P' X - (63) 
The background momentum-conservation equation (|^) gives the evolution equation of the field: 

^{r 2 V]A )^ A = ^ A A + 2v A ^ A = 0. (64) 
It is useful to notice that for scalar field matter 

t t r \A} A _2 A 2mfjawking , a r\ 

Vq = — = r - v A v = - -. (65) 

The scalar field has a perturbation J^i m YAip. We can construct a gauge-invariant perturbation as 

* = Atp-p c cp lc , (66) 
in terms of which the gauge-invariant perturbations of the stress-energy tensor are 

Tab = <&\ A ip\B + ®\b¥\a - 9ab(P\f^ F + ^9ABk DF <p\ D <fi\ F - \^.abV\f<^ ¥ \ (67) 

T A = $<p\A, (68) 

T 3 = kQ-<p lD & D + h DF <p lD <p lF , (69) 

T 2 = 0, (70) 

L A = 0, (71) 

L = 0. (72) 

Notice that there are no odd perturbations. 

Again, the momentum-conservation equation ( |3l"| ) gives the evolution equation for the matter perturbation, that is, 
the perturbed scalar wave equation: 

I > : ^(rH lA )\ A ^±11$ = i^fe^A)!* _ {k + v)lBip \B, (?3) 



Equations ( p2| ) and ( |33| ) are redundant for scalar field matter. If <p = 0, matter and metric perturbations decouple. 
To enforce regularity at the origin, we define 

$ = r'$, (74) 

where 3> is O(l) at r = 0. 



V. CHOICE OF FRAME AND COORDINATE SYSTEM 



In the remainder of the paper we shall use the radial basis defined by 

A 

n A = — W = v, U = 0. (75) 
v 

There is a system of coordinates naturally associated with this basis, which uses r as a coordinate: the familiar 
"Schwarzschild-like" coordinate system, in which the metric is 

ds 2 = -a 2 (t, r)dt 2 + a 2 (t, r)dr 2 + r 2 dQ 2 , (76) 

In these coordinates, the derivatives in the va frame take the form 
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-10/ 



"J fl -I" J 



at' 



<9r 



(77) 



This is not yet the coordinate system we shall use, but it is useful as an intermediate step in the presentation of the 
final coordinates. 

The Choptuik critical solution is a solution of the Einstein-real massless scalar field system defined by its self- 
similarity together with regularity. We introduce coordinates x and r adapted to self-similarity of the spacetime. The 
background solution has the geometric property of being discretely self-similar (DSS), which in our coordinates means 
that ip(r + A, x) — (p(T,x). The metric coefficients a and a defined in ( f76| ) have the same periodicity. Coordinates 
with this property are not unique. We make the following choice (in terms of the Schwarzschild-like coordinates): 



m [ -A 

ro 



(-1 



(78) 



where ro is an arbitrary scale. In the following we set it equal to 1. Our choice has the following properties. Surfaces 
of constant r coincide with those of constant t, and r increases with t. Therefore r is a good time coordinate, as well 
as being the logarithm of overall spacetime scale. The origin r = coincides with x = 0. We choose the function £o( T ) 
such that the past light cone of the point r = t = coincides with the surface x = 1. The domain of dependence of 
the disk < x < 1 on any spacelike surface is therefore given by < x < 1. We can therefore work on the numerical 
domain 0<a;<l,0<T<oo without requiring boundary data on x = 1. If we extended our perturbation initial 
data to x > 1, that part of the data could not influence x < 1. Wc can therefore determine exponential growth or 
decay on the domain < x < 1 alone. 

In these coordinates the frame derivatives in the radial frame are 



f = a- 1 e T 



df 
dr 



dr J dx 



J dx 



and the spacetime metric in these coordinates, but expressed through a and a, is 



ds 2 



l dr 



d£o 
dr 



■ dr 



-x 2 e 2 ^ dn 2 



(79) 
(80) 

(81) 



The background Einstein equations and a few more definitions are given in appendix 



VI. ODD PERTURBATIONS OF THE CHOPTUIK SPACETIME 



As we have seen, both La and L vanish and therefore the odd metric perturbations decouple from the matter 
perturbations. This implies [see equation (|6^)] that for I = 1 we have T = 0, and hence II = 0, if we demand 
regularity at the center. kA is then pure gauge. All / = 1 odd perturbations are therefore pure gauge. For I > 2 
equation ( |3l| ) is, in the radial basis, 



I > 2 : -n + n" + [v + 2(1 + l)v]W - - (I 2 - 4)V U - 
This equation is equivalent to the first-order system 



0. 



(82) 



(n)-+(n')' + 5i = o, 

-(u')-+(uy + s 2 = o, 

-(n)- + n = o, 
-(ny + fr = o, 



(83) 

(84) 

(85) 
(86) 



where 



Si = -/ill + \v + 2(1 + l)v}W - (I 2 - 4)Vbn, 

s* 2 = -/in' + vti. 



(87) 
(88) 
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Note that equations (J84[-p6[) are really identities that need to be added to the system when we consider it, H and IT 
as independent variables. From this first-order point of view, we now have three evolution equations (which contain 
dot-derivatives) and one constraint (which does not). The three characteristics are the light rays and the lines of 
constant r. Note that this causal structure is independent of any particular choice of coordinates. Now we introduce 
coordinates (t, x) . We also rescale II and its derivatives so that the rescaled variables u are precisely periodic in r if 
(and only if) the perturbed solution is DSS. Consider a perturbation of a self-similar background so that the sum of 
background and perturbations is again self-similar (to linear order in the perturbations). To find the scaling behavior 
of IT, we note that the tensor kAS a must scale like the metric itself. S a scales trivially, so that &a scales like the 
metric itself. On the other hand e AB scales like the inverse metric, so that the scalar e k^B scales trivially, that is, 
it is periodic in r for a DSS perturbation. Therefore n = r 2 ~ l e AB (r~ 2 kA)\B scales like r~ l , that is like e lT . We also 
note that each frame derivative adds a power e T . In order to cancel this scaling behavior, we define 

Ul = e - (z+1)r n, (89) 

u 2 = e- ( '+ 1)T II', (90) 
u 3 = e- lT IL. (91) 



The final form of the equations is then 



du . du „ . . 



where the 3x3 matrix A3 is 



with 



.l, = dia.i.Ll 2 .A„) i I ^ ) W) 



A2 = [-(°/a)e-*> Ao J' (94) 

We first consider the transport part of the equations. The characteristic speeds, or eigenvalues of A3 are 



a 



A , A± = A ± -e _i;o . (95) 
a 

Ao and A+ are always positive, while £0 has been chosen so that A_ changes sign at x = 1 by definition. That is, £o( T ) 
is defined by the equation 

1 -£)«--(:).-.■ (96) 

This definition means that for < x < 1 the characteristic speeds Ao and A + are positive, and A_ is negative. At 
x = 1, Ao and A+ are still positive, and A_ is zero. Therefore no boundary condition is required at the boundary 
x = 1, because no information crosses it from the right. At x = 0, all u are either even or odd in x, so that boundary 
conditions are obtained trivially. 

The source terms in the final equations are 

si = -ae- {l+2 ^ T Si + + l)ui = -a [-/mi + Du 2 + 2(1 + l)vu 2 - (I 2 - A)V u 3 ] + (I + l)u u (97) 

s 2 = -ae- (l+2)T S 2 + + l)u 2 = -a [-/2u 2 + vu Y \ + (I + l)u 2 , (98) 
S3 = — au\ + IU3. (99) 

We have used rescaled background coefficients that are periodic in r on the DSS background. Using the background 
Einstein equations they are 

u ee e- T n = -^-XY, (100) 
/ X 2 + Y 2 + 1 —f^), (101) 



e£°x 
1 

et°ax ' 
1 - a- 2 



(102) 

% = e- 2T V Q = t-^—. (103) 
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Note that at x = 0, these are regular except for v <~ x 1 . (The background fields X and Y are defined in the 
appendix.) Finally, the constraint equation becomes 

7jJ = ae ^- (104) 

As free initial data we can take u\ and U3, and we obtain u 2 by taking the derivative of 113. Numerically it is more 
stable to take u\ and u 2 , plus the value of u 3 at x = 0, as free data, and solve for u 3 by integration. 

VII. EVEN PERTURBATIONS OF THE CHOPTUIK SPACETIME 

A. General case I > 2 

The even perturbation equations are far more complicated. We discuss the cases I > 2, I = and I = 1 separately 
beginning with the general case I > 2. 

The vanishing of the matter perturbation T 2 makes fc^B traceless (r) = 0). Therefore the even perturbations are 
described by x, jp, k and <S. These obey the following set of equations: 



- I + $" - /xl + (v + 2(1 + !)«)$' + ^-/V + ^(F 2 - A 2 )^) $ 



-2(fc + r 2 x)i(F - ^ A) - 2?/>(A + (u - i/)Y) = 0, (105) 
-S + k" -n% + (v + 2(1 + l)v)k' - l 2 V Q k - 2x + 2 (v + ^(A 2 + Y 2 )^j (k + r 2 X ) + ^XY$ - ^A$ = 0, (106) 

- 8 

r 2 x + 2k + rip' + (2i/ + (Z + l)v)rip + 2^(k + r 2 x) + -F$ = 0, (107) 

- - - 8 - 

rip + r 2 x' + (l + 2)r 2 v X + 2v(k + r 2 x) + 2\irip - -A$ = 0, (108) 

+ (lV + ^Y 2 ^j (k + r 2 x) + ^(Y$ + X(lv$ + $')) + np + ^Y^j = 0, (109) 

-{%)' -(1 + 2)vk + iik' + 2)nvk - rvft - r^> {(I + l)v 2 +V + + -^(A 2 - Y 2 )^j - 2 l ivr 2 x 

--(X& + Y&) - -vY$(l + 2) = 0. (110) 

r r 

Again, we rescale the variables so that they are periodic in r if and only if the perturbed spacetime is still DSS: 

ui = e-^ T k>,u 2 = e-< I+1 > T $',u 3 = e~' T $, 
u A = e-^ T k,u 5 = e-^ T k',u 6 ee e- lT k, 

u 7 ee e-V+^rx, u s ee -e-V+^ii. (Ill) 
There are 8 evolution equations of the form 

du A du , „. 

d-r +A »0-x +S = °> (112) 

where the matrix As is 

A 8 = diag (A 3 , A 3 ,A 2 ), (113) 

and s is the vector 
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Si = (I + l)ui — a 



-Jim + + 2(1 + l)v)u 2 + -l 2 V + 



x 2 e 2£„ 



(Y 2 -X 2 ))u 3 



-2(e~ T Y -vX) 



1 



xe 



Co 



u 6 +u 7 )+ 2(e T X + (v - v)Y)u & 



52 = (I + 1)«2 - a[vui - /XU2], 

53 = Z«3 - aui, 
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S4 = (I + 1)«4 — a 



2 



«I« 3 - M«4 + (v + 2(1 + l)v)u 5 + (2 - l 2 )V + 



a; 2 e 2 £° 



(r 2 + x 2 )U 



+ ^(-a~ 2 + 2(X 2 + Y 2 ))u 7 - -j2u s 



S5 = (Z + l)t*5 — a[^U4 - ytm 5 ] 
s@ = Zwg — aii4, 

8 



S7 = (Z + 1)«7 — a 
s 8 = (Z + l)u 8 - a 
There are four constraints 



c 2 e 2£ 



x 2 e 2 £° 



Yu 3 



Xu 3 + 



2 



xe 



Co 



(w, 4 + [iu e ) - 2[iu 7 + (2f + (I + l)u)tig 
P«6 + (2P + (Z + l)v)u 7 — 2/2us 





du 3 
dx 


= ae^°u 2 




du e 
dx 


= ae^'us 


9a; 


b 7 
H u 7 

X 


= c 7 , 


du$ 
dx 


b% 

H «8 


= eg, 



where 
fe 7 = 



2 + 1 + 1 2 Z + l 



4F' 



&8 = a 2 
c 7 = a 
c 8 =a 



9^5 
9l*4 

9a; 



- a 2 e«° <! - 



a 2 e 6 ° 



a;e?° 
4 



(yui + X(u 2 + lvu 3 )) + Jim - 2(1 + l)vu 5 + l 2 V a + Ivv - 



x 2 e 2 t° 



M-6 



(Xui + Y(u 2 + (Z + 2)tw 3 )) — (Z + 2)wm 4 + fj,U5 + (I — 2)[ivu 6 [iu 7 

xe^o a 



(114) 

(115) 
(116) 



(117) 

(118) 
(119) 

(120) 
(121) 



(122) 
(123) 
(124) 
(125) 

(126) 
(127) 
(128) 
(129) 



The causal structure of the equations is similar to the odd case, because A% is constructed from A 2 and A 3 . The 
characteristics of A 2 are just the ingoing and outgoing radial null geodesies. u\, u 2 and u 3 on the one hand, and U4, 
u 5 and u 6 each form a wave equation with a mass-like term, while u 7 and us form a massless wave equation. The first 
two constraints are also identical to the odd perturbation case, and can be solved for u 3 and u§ by integration, or for 
u 2 and U5 by differentiation. Again we choose the former in the numerical treatment, taking the value of u 3 and uq 
at x — as free initial data, together with u\, u 2 , u 4 and u 5 . 

The next constraint equation contains u 7 but not u 8 , and is therefore a linear ordinary differential equation (ODE) 
for u 7 . Once u 7 is known, the last constraint can be solved as an ODE for u s . We solve these ODEs by a second-order 
implicit method, in order to finite-difference all constraints in the same way. Both the evolution equation for u 7 and 
the constraint for u s require the following condition at the origin x — for all r in order to be consistent: 



2 Ui + 4-?" u 3 - (Z + l)u 8 = 0(a; 2 ). 
e«° 9a; 



(130) 



We solve this constraint for the value of u$ at x = 0. The value of u 7 at x = is zero by definition. These boundary 
conditions complete the constraints for u 7 and u$, which are then determined completely, given u\ to u§. 
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B. Special case / = 



For I = a general perturbation is described by the objects (Jzab > k, Tab, T 3 ), which are not gauge-invariant: under 
an arbitrary coordinate transformation generated by the vector ^dx^ = l;AYdx A these objects change as 

kAB — ► k A B - (£a|B + £b\a), (131) 
k — ► fc - 2v A i A , (132) 
Tab — > Tab ~ ^ab\c\ ~ tAci°\B ~ *bc£ |a> (133) 

T 3 ^T 3 -l(r 2 Q)| D e D - (134) 

Therefore we have to impose two gauge conditions. In our case we want to maintain the form ( |76| ) of the metric 
during perturbation, so we perform a gauge transformation to obtain fc = tp = 0. Then, metric perturbations are 
described by ij and \. By regularity they are 0(1) and 0(r 2 ) at the center, respectively. The condition k — fixes the 
projection of £ on v A completely, but -0 = fixes the orthogonal part only up to a residual gauge freedom ^a — Jua 
where the scalar / obeys the equation /' = vf . This latter equation can be thought of as an ODE in r at constant t. 
We can give the boundary value for this ODE at each moment of time, so the residual gauge is an arbitrary function 
of time. We use it to set rj = at the center. 
Using @, @ we define 

r) = fj, x = (f) + r] = r 2 x, (135) 

where \ IS 0(1) at the center, but fj is 0(x 2 ), due to our gauge choice. The scalar field perturbation <I> is already 
0(1) and even at the center, compare J74|). Equations (|45|-[47|) and ( f73| ) are then 

-X' + (1 + 2a- 2 ) X + 4F 2 ( 4 - x) ~ + X &) = °> (136) 

a \r z / r 

— fj' + 4Y 2 (^-x) --{Y& + X&) = 0, (137) 
ar \r z / r 

-%--(X$ + Y&)=Q, (138) 
a r 

-$ + $" - -oir$ + — (1 + 3a- 2 + 2X 2 - 6Y 2 W 
r 2r 

+ -fj + aX(l - 4Y 2 )x + ^aXY 2 fj + 2 ( % - X ) rY = 0. (139) 

The last equation is the wave equation for the matter perturb ation . We do not have an evolution equation for fj. 
Instead, we ha ve to calculate it by integration of th e co nstraint (137). Finally X can be calculated from the evolution 
equation (138) or by integration of the constraint (136). 

Again we rescale the variables. We also reorganize the variables to eliminate fj from the equations. (This is the 
same trick as using Y instead of tp to simplify the background equations.) 

it! = e~ r - ^T] \ , u 2 = e _r $', u 3 = $, 

W4 = fj,U5 = e~ T rx. (140) 

Variables 112, M3, u 5 ) verify the following evolution equations: 

du A du , , . 

— +A 4 — + s = 0, 141 

OT OX 

where the matrix A4 is 

A 4 = diag(A 3 , A ), (142) 

and s is the vector 
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Si = ul — a 



Go 



XYu x + -v 



a ,1 



(-+X 2 -3Y 2 ) ua + (-2a 



s 2 = u 2 - a 



aX 

xe^° 
«o _ 



f 1 - AY 2 



xe^o 2 
2e~ T Y } u 5 



XY 2 

c 2 e 2{ 



3 T -^)w 4 

xeM> 



S3 



2 
Y 



Ml 



S5 = Ub - OL 



xe^° 
4a 



-Vb 



«4 



(X 2 + 5Y 2 ] 



2aXY 



Ul 



-U2 



2a- 



X 2 Y 



xe€° 

There are three constraints: 



XY 

(Xui + Yu 2 ) + 4a , 0Cn m 



x 2 e 2 Z° 



du 3 
dx 

dx 



ae^°u 2 . 



c 4 , 



where 



du 5 b 5 

— 1 u 5 = c 5 , 

ox x 



c 4 = 4a 2 e« (Fui + Xu 2 + Y 2 u 5 ) 
b 5 = l + a 2 (l-4r 2 ), 
4a 2 

c 5 = (Yui + Xu 2 ). 



X 



U4 



4aF 3 

xeA° 



(143) 

(144) 

(145) 
(146) 

(147) 
(148) 
(149) 



(150) 
(151) 

(152) 



Note that we have a constraint, but no evolution equation, for 114. We have in fact constraints for U3 and both U4 
and U5, so that the only degrees of freedom are those of a wave equation. We obtain 114 by solving a constraint at 
each time step, starting from the gauge condition U4 = at x — 0. 



C. Special case / = 1 

For / — 1 a general even perturbation is described by the objects (kAB, k, Tab, Ta, T 3 ), which are only partially 
gauge-invariant: under an arbitrary coordinate transformation generated by the vector ^dx^ — ^AYdx A + r 2 £Y. a dx a 
these objects change as: 

k AB — » k AB + (r 2 ^\A)\B + {r 2 i\ B )\A, (153) 
k^k + 2Z+(r 2 )\% A , (154) 

Tab — > Tab + r 2 (t AB \c^ C + t AC S} C B + t B c( lC A ), (155) 
T A ^T A + r 2 {t AB S} B - Q£\ A ), (156) 
T 3 — ► T 3 + 2Q£ + {r 2 Q)\ A t;\A- (157) 

We see that there is invariance under the £a part of the transformation. Therefore we have to impose just one partial 
gauge condition. The most interesting gauge condition is k = 0, because then we can eliminate all second derivatives 
from equations (|45| - pl|) . Now matter perturbations are described by T),ij),X, which are 0(r), 0(r 2 ) and 0(r 3 ) at the 
center, respectively. The condition k = does not fix the gauge completely, and again we have a residual gauge 
freedom of functions £ obeying equation r 2 vt;' + £ = 0. We use this freedom to set 77 ~ 0(r 3 ) at the center. 
Using @, (§§, (§|) and @ we define 

rj = rfj, -0 = r 2 V>, \ — 4> + V — r ^X, ^ = (158) 

where the barred variables are even and 0(1) at the center, except fj, which is 0(r 2 ), due to our gauge choice. 
Equations (Hffc]) and @ are then 



13 



V 4 fX 



-$ + y$ + X$' =0, (159) 



-X' + (2 + 3a- 2 -AY 2 )x + W 
—f)' + 2(1 - a' 2 +X 2 + Y 2 )% + (a- 2 - 2X 2 - 2Y 2 ) X - -XY$ + = 0, 



ar 



(160) 



-J) 1 + (2 + a- 2 + 2X 2 - 2Y 2 )tP + ArXY (x ~ ^) + 4 ^$ + X$ + F$' j = 0, (161) 

J, 4 / Y - - - \ 
' — \ 1 i '->• ,;, - /!•,' 1 — i , (162) 

(163) 



-% - (1 - 4F 2 )- - - ( — $ + x$ + ) = 0, 

a r r \ar 



-V - 4XYip + (3(1 - a" 2 ) + 2(X 2 - Y 2 ))^ - (1 - a" 2 + 2(X 2 - Y 2 ))r X + 4 ( — $ + y# + ) = 



3X 



ar 



') 



-$ + $" _ -aXY® + (1 + 7a- 2 + 2X 2 + 2F 2 ) — + -V + 8— $ 



y 1 



+-»? + (aX - 2rY)x + {aX + 2rY)% 



aY 



2r 



(1 - 3a~ 2 - 2X 2 + 2Y 2 ) 



- 2X) tp = 



tp = 0. (164) 



Again we rescale and regroup the variables: 



Ul = e- 2T 



Y 



<f> - — rj I , U2 = e 2r $', u 3 = e T $ 



-2t„ 



?i4 = e 'r],U5 = e ~'rx,UQ = ~e T ip. 
The variables (ui, u 2 , u 3 , «5, u 6 ) obey the following evolution equations: 



where the matrix A§ is 

and s is the vector 

2aXY 



du , du 

-7T + A b — + S = Q, 
OT OX 



M = diag(A 3 , A ,A ), 



(165) 



(166) 



(167) 



s\ = 2ul — a 



aX „ _ • \ / , aF 



8F 2 



M4 



s 2 = 2u 2 - a 



2e- T Y )u 5 +[3vY r (l - 2X 2 + 2Y 2 ) + 2f- T X u 6 



-V0 + -J- (X 2 + Y 2 ) Ul 



2aXY 



(168) 



xe 



-u 2 



8XY 

x 2 e 2 t° 



u 3 + -2aY(V Q + 



Y 2 



x 2 e 2 Z° 



) + e" 



xe 



So 



1*4 



+ - 



Y 



xe 



Co 



-- + 2a(X 2 + Y 2 ) )u 5 - 



AaXY 2 



xe 



to 



«3 = "3 - 

s 5 = 2u 5 - a 
sq = 2uq — a 



Y 

ui H ^it4 

SO 



xe 



4a (Xui + y« 2 ) + zJ^Ym + 4a^^ U4 - ^(1 - 4F 2 ) U6 



xet° 
4a 



(Ym +Xu 2 ) 



x 2 e 2 £° 
12 

c 2 e 2Co 



Xu 3 + (3aVb + 



2a 



x 2 e 2£ 



(X 2 + Y 2 )) Ui 



- ( axe^Vo + ^ r {X 2 -Y 2 )) u 5 + ^XYu 6 



2a 



xei° 



4a 



xei° 



(169) 
(170) 
(171) 

(172) 



There are four constraints: 



du 



dx 



- = ae^°u 2 , 



(173) 
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8u4 &4 

-5 1 u 4 = c 4 , 

ax x 


(174) 


du 5 b 5 

— 1 u 5 = c 5 , 

ox x 


(175) 


du 6 b 6 

— 1 u 6 = c 6 , 

Ox x 


(176) 



where 

b i = -2 + 2a 2 (l+X 2 + Y 2 ), (177) 
& 5 = 2 + 2a 2 (l - 2F 2 ), (178) 
b 6 = l + 2a 2 {l + X 2 -Y 2 ), (179) 
8aA 

c 4 = u 3 + e ?0 (-l + 2a 2 (A 2 + Y 2 ))u 5 - 4e^a 2 XYu 6 , (180) 

4a 2 

c 5 = (Y Ul +Xu 2 + vXu 3 ), (181) 

4a 2 

c 6 = (X Ul + Yu 2 + 3vYu 3 + XYu 5 ). (182) 

a; 

Note that again we do not have an evolution equation for 114, and that we have constraints for all variables other than 
Mi and u 2 , so that the only degrees of freedom are those of a wave equation. There is a consistency condition at the 
center: 

4 dY , 2x , N 

u e - ——u 3 = 0(x 2 ). 183 
e«° ax 



Again we impose U4 — at x = as a gauge condition. 



VIII. NUMERICAL RESULTS 



Our numerical code, and the tests we have performed, are described in an appendix. Here we only summarize three 
important points. 

The code treats the boundaries x = (center of spherical symmetry) and x — 1 (boundary of domain of dependence) 
in exactly the same way as all other points. On a flat empty background spacetime, it is second-order convergent, the 
origin x = is stable, and waves cleanly leave the computational domain at x = 1 without numerical backscatter. 

On the Choptuik background we observe second-order convergence for most values of x and r. Convergence of 
a lower than second order is observed near x = 0, twice per period in r. These are the values of r where certain 
coefficients of the background solution change rapidly in time, namely at the minima and maxima of the background 
scalar field. A typical solution (as we shall discuss below) is an exponentially damped quasiperiodic oscillation. 
Convergence inevitably breaks down at large r for two reasons: the oscillations at different numerical resolutions 
gradually drift out of phase, and small differences in the exponential decay rates at different resolutions have a 
cumulative effect on the amplitude. 

As we discuss in detail in the appendix, the numerical code has a subtle instability which becomes apparent only at 
high I at high resolutions. The instability is already present in the free wave equation (in self-similar coordinates) on 
Minkowski space. We have found a way of repairing it in Minkowski space, but it persists on the Choptuik background. 
At low resolution this instability can be neglected, and we see convergence up to a resolution of Aa; = 1/800. 

In spite of the inevitable absence of pointwise convergence at late times, and in spite of the numerical instability, 
our main result appears secure: all non-spherical physical perturbation modes, for all initial data, decay exponentially 
in t. The exponential decay is typically rapid. Only for even I = 2 perturbations is the decay quite slow, but there 
(as for low I in general) we have good convergence of the solution itself, and therefore the decay exponent. 

Due to the discrete self-similarity of the background solution, the perturbations decay in a complicated fashion, with 
the exponential decay apparent only over many periods. The background-dependent coefficients of the perturbation 
equations are periodic in t (at constant x). Therefore the general form of the perturbation is a sum of terms of the 
form 

u(x, t) = Re [C e XT f(x, r)] , (184) 
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with C, A and f(x,r) all complex, and f(x,T + A) = f(x,r). Once the most slowly decaying mode dominates, only 
one such term is left. In real notation, it is 



u(x,t) = e KT [C 1 cos( W r)/ 1 ( a :,T) + C 2 sin( W T)/ 2 (x,r)], (185) 

with re = ReA, to = ImA, C%, C%, fi and f% now real, and f\ and f% again periodic. This means that u(x, r), even after 
the exponential decay has been taken out, is not periodic in r unless u> is commensurate with 2n/ A. Furthermore, 
C\ and C2, and in particular their ratio, depend on the perturbation initial data. Therefore, the complex exponent 



A is not easy to read off. Nevertheless, to the extent to which they are approximated by (185), the Fourier transform 



in r of the data with the exponential decay taken out should be peaked around the set of frequencies 

N—+LU (186) 

for integer N. The background is not only periodic in r with period A, but has an additional symmetry. The 
background scalar field obeys </?(t + A/2, x) — —(p(r, x), while the background metric coefficients obey g{r + A/2, x) = 



g(r,x). The perturbations inherit this additional symmetry. Therefore, in the spectrum (186) of the scalar field 
perturbations u\ to 113, only odd integers N appear, while in the spectrum of the metric perturbations 1/4 to u&, only 
even value of N appear. This must be taken into account when we read off u> from the spectrum. Because the N 
are either even or odd, loA/2tt is defined modulo 2 (and not modulo 1 as one might expect), and we define it to be 
< uiA/2ir < 2. For example, with the highest peak in the spectrum of Ui at loA/2tt — 5.3, and the highest peak in 
the spectrum of M4 at luA/2tt = 6.3, we consistently obtain luA/2tt = 0.3. 

The only exception to this complicated behavior are the spherical (I — 0) perturbations. At large r they are 
dominated by a single growing mode with A real. (The fact that there is a single growing mode is of course at the 
center of critical phenomena in critical collapse, and this unique A must then be real because the background is real.) 
Here we can read off both re and f(x, t) quite clearly. We find kA — 9.21. This corresponds to re = 2.67, and a critical 
exponent for the black hole mass of 7 = 1/re — 0.374. This agrees to all three digits with the value of the critical 
exponents obtained from collapse simulations 0, and a perturbative calculation JT(| that is completely independent 
from the present one. 

For I > 0, we have obtained estimates of re and to by first adjusting the value of re until the rescaled lt re scaled = £~ KT u 
appeared to be neither increasing nor decreasing over many periods. The resulting u reS caicd is then quasiperiodic. We 
have carried out a discrete Fourier transform on the time series wi, rescaled (0, t) over a range of 10A. The result has 
sharp peaks spaced at intervals 47r/A due to the additional symmetry in the background mentioned above. In the 
special cases I — and I — 1 the function u\, rescaled is clearly periodic (lu — 0), and the line spectrum is very sharp. 

The estimated values of re and lu are tabulated for different resolutions in Table 1. As an example, we show the 
value of Mi for even I = 2 perturbations at x — as a function of r, after an exponential decay has been taken out, in 
Fig. [|. In Fig. || we show the low frequency part of the discrete Fourie r tra nsform of Fig. [|. The quasiperiodic nature 



of the signal becomes clear in that there is a series of peaks obeying (186). 

As the background spacetime is periodic in r, and the perturbation equations are linear, evolving the perturbations 
for one period is equivalent to multiplying them by a transfer matrix. For odd perturbations, this matrix has size 
(2N) 2 , and for even perturbations (4iV) 2 , where N is the number of grid points in x, and two and four respectively is 
the number of degrees of freedom. For N = 50 and 100, we have verified that the logarithm of the largest eigenvalue 
of the transfer matrix agrees with AA. These matrices contain of course all the information that there is about the 
system, but for larger N the computation time and memory requirement for calculating these matrices and their 
eigenvalues quickly becomes prohibitive, scaling as A 4 . However, if we use generic initial data, in which no u vanishes 
at any x (except odd u at x = 0), we have a mixture of all perturbation modes, and at late enough times the most 
slowly decaying mode has taken over. 

With increasing I, the even parity numerical code appears to be more and more sensitive to small errors in the 
background solution, as the solutions obtained at different resolutions drift apart more and more rapidly. The solution 
at late times depends sensitively on the initial data, so that the system looks chaotic. This problem may be unavoidable 
with any code. Our results still seem to capture the correct overall behavior, as the values of re and to obtained at 
different resolutions differ much less than the actual time series. We believe that the explanation is that different 
resolutions agree reasonably well on the periodic functions f\ and f2, but that the initial data-dependent coefficients 
C\ and C2 take essentially random values at late times for different resolutions. 

In summary, we find that both even and odd perturbations decay exponentially for all physical values of I. It is clear 
that perturbations with large I will decay more and more quickly because of the presence of the terms u iT = —lu + . . . 
in all evolution equations. (These terms are introduced by the scaling of perturbations with r l to keep them regular at 
r = 0.) The numerical evolutions confirm that higher / modes decay more and more rapidly. We can therefore affirm 
that all values of I decay, even though we have checked this explicitly only for the lowest few values. The most slowly 
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01 23456789 10 

FIG. 1. ui versus r at x — 0. An overall exponential decay has been compensated for. The scale on the vertical axis is 
irrelevant, as the equations are linear. On the horizontal axis we have marked background periods, that is r/A. The sharp 
peaks are typical features. Although it is not clear from this plot, they are perfectly smooth. 
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2 4 6 8 10 12 14 16 18 20 



FIG. 2. The low-frequency end of the discrete Fourier transform of the previous figure. The vertical scale is again irrelevant. 
On the horizontal axis we have marked frequency in units of the background frequency, that is (a>A)/(27r). The quasiperiodic 
nature of the signal shows up in the peaks situated at (ojA)/(2tt) = 1.3, 3.3, 5.3, . . .. The spectrum decays rapidly at high 
frequencies. 
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System 






(kA,cjA/2tt) 






grid points 


100 


200 


400 


800 


1600 


even I = 


9.24, 0.0 


9.21, 0.0 


9.21, 0.0 


9.21, 0.0 


9.21, 0.0 


even I = 1 


-0.34, 0.0 


-0.31, 0.0 


-0.48, 0.0 


noisy 


-0.30, 0.0 


even I = 2 


-0.08, 0.3 


-0.07, 0.3 


-0.06, 0.3 


-0.07, 0.3 


-0.07 0.3 


even I = 3 


-1.63, 1.6 


-1.65, 1.6 


-1.65, 1.6 


-1.65, 1.6 


-1.66, 1.6 


even I = 4 


-2.8, 0.9 


-2.9, 0.9 


-2.9, 0.9 


-3.0, 0.9 


noisy 


even 1 = 5 


-4.0, 0.2 


-4.25, 0.2 


-3.9, 0.2 


-3.65, 0.3 


noisy 


odd Z = 2 


-2.20, 1.9 


-2.28, 1.9 


-2.30, 1.9 


-2.30, 1.9 


-1.8, 3.0 


odd I = 3 


-3.13, 1.3 


-3.23, 1.4 


-3.27, 1.4 


-3.28, 1.4 


noisy 


odd I = 4 


-4.05, 0.7 


-4.20, 0.7 


-4.25, 0.7 


-4.27, 0.7 


noisy 


odd Z = 5 


-5.0, 0.0 


-5.2, 0.0 


-5.2, 0.1 


-5.3, 0.1 


-5.3, 0.1 



TABLE I. Summary of eigenvalues A, read off from u\(x = 0,r). The values of kappa were obtained by eliminating an 
exponential factor. The values of ui were obtained from a discrete Fourier transform of the result. Values of uo/S./2-n are defined 
modulo 2, while peaks in the discrete Fourier transform of ui are located at ll>A/2tv + 1 + 2N for non-negative integer N. As 
we have integrated over a range of 10A in r, lo/S./2-k can only be estimated in multiples of 0.1. Results marked "noisy" are 
dominated by numerical error. In the Fourier transform this shows up as high frequency noise. 



System 






(kA,cjA/2tt) 






grid points 


100 


200 


400 


800 


1600 


even I = 2 


1.1, 0.3 


0.15, 0.3 


0.0, 0.3 


-0.05, 0.3 


-0.07, 0.3 


even I = 3 


0.25, 1.4 


-1.3, 1.6 


-1.55, 1.6 


-1.63, 1.6 


-1.65, 1.6 


even I = 4 


-0.53, 0.5 


-2.0, 0.8 


-2.75, 0.9 


-2.8, 1.0 


noisy 


even I = 5 


-1.45, 1.5 


-2.77, 0.0 


-3.3, 0.1 


-3.2, 0.2 


noisy 


odd I = 2 


-1.74, 1.8 


-2.18, 1.9 


-2.30, 1.9 


-2.30, 1.9 


noisy 


odd I = 3 


-1.95, 1.2 


-2.95, 1.3 


-3.2, 1.4 


-3.25, 1.4 


noisy 


odd I = 4 


-2.2, 0.4 


-3.55, 0.7 


-4.1, 0.7 


-4.25, 0.7 


noisy 


odd I = 5 


-2.7, 1.4 


-4.1, 1.9 


-5.0, 0.0 


-5.2, 0.1 


noisy 



TABLE II. The same with an alternative finite differencing method for I > 2, defined by Eq. (B9). Note that convergence 
is much slower, but that for 800 grid points these results agree quite well with those of the other method. 



decaying mode occurs in the I = 2 polar perturbations, with A ~ —0.07 x (1/A) + 0.3 x (2iri/A) ~ —0.02 + 0.55i. 
As this mode decays so very slowly, there may be an intermediate range of p — p* for a given one-parameter family 
of initial data where this perturbation becomes universally visible. For p — p* small enough, however, the spherical 
universal solution will again dominate. 
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APPENDIX A: BACKGROUND SOLUTION 



Following Choptuik, we introduce the scale- invariant and dimensionless auxiliary fields for the spherically symmetric 
scalar field system as 

X = V2^rtp', Y = V2^np, (Al) 

in the radial basis. (Do not confuse Y with a spherical harmonic.) We use the dependent variables g — a/a, a, X 
and Y, and the coordinates x and r. With the shorthand 

D = xge^{l-i Q , T ), (A2) 
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the wave equation for (p is equivalent to the first-order system 



X.. 
Y, 



The Einstein equations are 



1 



1 D 



1 — D 2 \ D 1 



-[±(1 + a 2 ) + a 2 {X 2 - Y 2 )]X + gxe^Y T 
[i(3 - a 2 ) + a 2 (X 2 - Y 2 )]Y + gxe^X, T 



.9(1 -a 2 ), 

-[l~a 2 + 2a 2 (X 2 + Y 2 )}, 



-(1 - £o,r)xa, x + 



2a J 



gxe 



XY. 



(A3) 

(A4) 
(A5) 

(A6) 



In order to exclude a conical singularity at r = 0, we impose a = 1 at x = 0. In order to fix the remaining coordinate 
ambiguity t — > fit) we impose g = 1 at x = Q. We make i = 1 an ingoing null surface by imposing 13 = 1 at x = 1. 
£o is not initially known, but is determined together with the dynamical fields X, Y, a and g of the critical solution. 

We have recalculated the background using the numerical code of Gundlach |l(J , slightly modified to use x instead 
of ( = lnx, which results in a better treatment of x — 0. If the solution is regular, X and Y vanish at x = 0. 



Therefore we work with X 2 



l X and Y 1 



l Y. 



and x = 1 are regular singular points of the equations. 



The regularity condition (vanishing of the numerator in the wave equation) is 



X 2 = \e^° [Y ltT + (1 - d$a/dT)Yi] 



at x = 0, while at x = 1 it is 



[1 + a 2 (l + 2Xl - 2^i 2 )] A 2 + [-3 



a 2 (l - 2Xl 



2y 1 2 )]y 1 -2(i-^| 



dYi 



dX 2 
~~dr~ 



= 0. 



(A7) 



(A8) 



The discrete self-similarity of the background is equivalent to periodicity of X, Y , a and g in t, w ith a period A 
that is initially unknown, a is treated as a functional of A, Y", 5 and £o 5 by solving equation (A6), with periodic 
boundary conditions in r for each value of x. Note that this equation is linear in a~ 2 . Periodicity is impo sed by 
representing X, Y, g and £0 through a (truncated) Fourier series, r-derivatives are calculated, and equation (A6) is 
solved, in Fourier space. This makes the numerical method a pseudo-spectral one. The y-derivatives are implemented 
through finite differencing on a grid equally spaced in x, and are solved by relaxation, together with the algebraic and 
ODE (pseudo-algebraic) boundary conditions at x = and x = 1. 

We have calculated the background solution using points 51, 101, 1601 on the range < x < 1, always with 128 
points per period < r < A. It was shown in |10| that this resolution in r is large enough so that numerical error is 
dominated by resolution in x and systematic error effects at x — and x = 1. 

We observe second-order convergence, measured by the maximal and root-mean-squared differences of X 2 , Y\, a and 
<7, from 51 to 1601 grid points in x (with 128 Fourier components in t). A and £0 (i n the maximum and root-mean- 
squared norms) also show convergence, but not to a distinct order. This is illustrated in Fig. [| For the perturbation 
calculations we have always used the high-resolution background, downsampled as necessary. 



APPENDIX B: PERTURBATIONS NUMERICAL METHOD 



The perturbation equations are of the form 

du ,du , . 

— +A— + Bu = 0, Bl 

OT OX 

where u is a vector of unknowns, and A and B are background-dependent matrices. 

As by definition x — 1 is an ingoing spherical null surface, the domain of dependence of perturbation data at r = 0, 
< x < 1 for r > is precisely < x < 1. Both scalar and gravitational waves travel only from smaller to larger x for 
x > 1. In order to implement a time evolution without an artificial boundary condition at x = 1, we use an evolution 
scheme that explicitly uses the characteristic speeds and therefore changes over smoothly to upwind x-derivatives 
for x > 1. The numerical method we have used is similar to that used for the perturbations of the perfect fluid 
critical solutions in a previous publication P]l3[|, but is second order in space and time. It uses second order one-sided 
derivatives in x, and is Runge-Kutta-like in r: 
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1e+01 f 




1e-04 ' 

100 1000 



FIG. 3. Convergence of the background solution with increasing resolution in x. On the horizontal axis we plot the logarithm 
of the number of grid points, on the vertical axis the logarithm of the difference between one numerical solution and the one 
with half the resolution. The two thick lines are the maximal and root-mean-squared differences (over both < x < 1 and one 
period in r) of all fields. The thin lines are the differences of the "eigenvalue" £o(t). The dashed line is the difference in the 
eigenvalue A. 
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n+i 



n+1 



At 
At 



4 <4 



2Aa 



3< 4< 



X j-2 



3u« 



2Ax 



B>1 



*j+2 



3u 



2Ax 



4u 



™+2 o r ' 



2Ax 



B 



n+i n+i 



(B2) 
(B3) 



Here A + + A- = A. In order to use the characteristic speeds in the finite differencing scheme, it is necessary to split 
A into a sum over its eigenvalues according to their sign, so that, for example, 



-43+ = 



' 1/2 -1/2 
-1/2 1/2 



'0 



+ A 



A a - = X- 



0, 



1, 




(B4) 



For x > 1, A_ becomes positive, so thatA + = A and A^ = 0, so that we do not need the downwind derivative there. 
We go just one grid point beyond x = 1, so that the last grid point just before x = 1 still has two points to its right 
in order to take a right derivative there. All grid points further to the right only require left derivatives. This means 
that we could have extended the numerical domain to large x. We have chosen < x < 1 + Ax because it is the 
smallest numerical domain in which we stay in the domain of dependence of the perturbation initial data for all r, 
while using a one-sided three-point stencil. (If we used a first-order differencing scheme, with two-point stencils, the 
numerical domain < x < 1 would be sufficient.) 

We might refer to the method just described as the second-order characteristic method. It is explicit and second 



order. One obtains an implicit method by averaging u n and u n+1 to obtain a new improved estimate for 



71+1/2 



-1/2 _ 



= {u n +u n+1 )/2, 



(B5) 



and iterating the pair (B5 B3) of equations until u n+1 has converged. Let us call this the iterated characteristic 
method. 

The boundary x — does not require special treatment, as u(—x) = ±u(x) for all u, so that ghost grid points with 
x < are available for taking derivatives. The one-sided differencing operators do not give exactly zero at x = even 
if analytically du/dx(0) = 0, but that is consistent: all terms in the finite difference equations combine so that u(0) 
remains zero to machine precision if u is odd initially. This also ensures that source terms of the form u/x for odd u 
in the evolution equations are well behaved numerically 

Another promising candidate for numerical treatment of the equations is the iterated Crank-Nicholson method. In 
this method, we need to treat the boundaries specially. At x — we have tried updating the boundary point by 
extrapolation from its next neighbors at each iteration, taking into account that u is either odd in x (and vanishes 
at x = 0) or even in x. Alternatively we have used the exact value du/dx — for even u and the finite difference 
stencil du/dx = \u{Ax) — u(0)]/Ax, which is second order at x — for odd u. At x = 1 + Ax, we have either used 
extrapolation, or the one-sided (left side only) second-order stencil of the characteristic methods. 

The constraints are solved by integrating from x — out to x — 1. For stability, we do not evolve any variable for 
which we have a constraint, but instead calculate it from the constraints, including at the half-step n + |. For simple 
integrations du/dx — c with u an even function of x, and u(0) given, we use the trapezoidal rule 



Aa; 

U j+ i = Uj + — (Cj + c i+i ) . 
For the ODE du/dx + bu/x — c, with b and c even in x and, therefore, u odd in x, we use 



u j+ i 



Ax bj+i 
Xj + Xj-\~\ 



1 - 



Ax bj 

Xj + Xj+l 



Ax 

— (Cj+C j+1 ) 



(B6) 



(B7) 



This scheme is second-order accurate at all grid points. For those variables u that are even in x and of 0(1) at x = 0, 
because c is odd and of 0(x _1 ), we use the same scheme, with coefficients 6—1 and cx, in order to first calculate the 
odd function u — ux. Then we divide by x to obtain u. In this case we extrapolate twice: first cx to x = 0, and then 
u/x — u to x = 0. 

The perturbations were calculated at different resolutions, related by grid refinements by a factor two in both x 
and r. As our lowest resolution, we used Ax = 0.01, with a Courant factor of At/Ax ~ 0.05. (The exact Courant 
factor is chosen so that the number of time steps for integrating the perturbations is a multiple of the number of 
time steps in the background, per period.) This small Courant factor is necessary for stability, apparently because 
some coefficients of the perturbation equations, although smooth, have very large gradients in r near x = 0. We have 



22 



also verified that the effect of using an even smaller Courant factor is negligible. Our highest resolution was finer 
by a factor of 16 in both space and time. The background coefficients a, a, X, Y were given in Fourier coefficients 
at 128 points per period in r and the required intermediate values of t were obtained by local cubic interpolation. 
We chose local cubic interpolation as it is much faster than Fourier interpolation, and because of limited computer 
memory. The interpolation is formally second order accurate, and all background coefficients are well resolved at 
this resolution. To separate the convergence of the perturbations from that of the background coefficients, we used a 
background solution obtained with 1601 points in x, and downsampled it by factors of 1, 2, 4, 8 and 16. 

As a test, we used all numerical methods on the trivial background of flat empty spacetime without a scalar field. On 
this background, the even matter and metric perturbations decouple. In fact, the even matter perturbation equation 
is identical to the odd master equation, and both are identical to the free wave equation. We are therefore testing 
the code on the free wave equation, with angular dependence Yi m , and in self-similarity coordinates, on the domain 
of dependence < x < 1. 

Already in this simple test, we find that the iterated Crank-Nicholson method with any of the boundary treatments 
discussed is unstable. The instability does not have a continuum limit in space. In fact, it changes sign about every 
grid point in space, and grows twice as fast in time when Ax is halved. Nevertheless, it appears to have a continuum 
limit in time. The instability changes smoothly from one time step to the next, and in fact, it is practically unchanged 
if At is reduced by a factor of 10 (at constant Ax). The instability is most apparent at x = 0, but in an implicit 
scheme, all grid points in space are of course linked. 

As we have not found a boundary treatment in which iterated Crank-Nicholson is stable, we now restrict ourselves 
to the two characteristic schemes. In flat space we find that they give essentially the same solution. Both are stable 
and second-order convergent for a long time. In particular, x — and x = 1 are perfectly normal points, as expected 
from the construction of the numerical scheme. At high resolutions, high I, and late times (for example, Ax = 1/1600, 
I = 6, t > 1) we see an oscillating instability in r near x = that leads to a breakdown of convergence near x = 0, and 
which blows up for sufficiently large / and high resolution. This instability appears to be a solution of the continuum 
equations, but for initial data which are provided by finite differencing error at late times. To demonstrate this, we 
have evolved a narrow Gaussian pulse originally centered around x = 0.5 in flat spacetime. After this pulse has left 
the computational domain, nothing should happen physically, and we are left with pure numerical error. We then 
extracted new Cauchy data at a late time, at a high resolution, and restarted these data with different resolutions 
(in space and time), down-sampling the error data as necessary. For some time we clearly see quadratic convergence, 
until new numerical error takes over. 

The instability is present already in the free wave equation in flat space, which in our rescaled variables and 
self-similarity coordinates is 

du x dui du 2 21 + 2 n 

-t— + x— u 2 -{1 + 1 Ml = 0, 

or ox ox x 

du 2 du 2 du\ . „ 

---o 1 ~ + 1 )«2 = 0. (B8) 
or Ox ox 

Recall that for any /, m is an even function of 0(1) of x and u 2 is an odd function of O(x). As the instability is 
centered at x = and becomes worse with increasing I, it must be linked to the term u 2 /x. We have generalized 
a well known trick for the spherical wave equation which consists in absorbing this term into the x-derivative. We 
rewrite the equations as 

We then finite difference u 2 always in the way suggested by the equations, using left and right second-order one-sided 
derivatives to obtain the characteristic method outlined above. Note that we only ever use the new derivative of 
U2, never the straight derivative du 2 /dx. The generalization to the even and odd parity perturbations of scalar 
field collapse, applied to the variables u 2 and u 5 , is straightforward. In particular, the coefficients of du 2 /dx and 
(21 + 2)u 2 /x, although now different from unity, remain equal to each other. All other variables u are differentiated 
directly with respect to x. When we use this finite differencing method for the flat space wave equation, the late-time 
solution that is pure numerical error is now smooth and decays exponentially instead of blowing up, as we had hoped. 
On the Choptuik background, however, the instability at the center is still not suppressed. Before the instability takes 
over, the two methods clearly converge towards each other. 

We have also tested convergence of the code on the (numerically generated) Choptuik background. Here, stability 
requires a smaller Courant factor. The differences between different resolutions are peaked at x = and are smooth 
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functions of x. For most values of r and x convergence is clearly second order, with the exception of certain value 
of t near x = 0, where the background coefficients are particularly rapidly varying. Here convergence still occurs 
(differences between resolutions decrease with resolution), but is not of a clear order: lower than second order for 
low resolutions and higher than second order for high resolutions. Somewhat surprisingly, second-order convergence 
disappears and then reappears many times. Apparently error is not just growing with time, but depends very strongly 
on the background. The explicit and iterated characteristic methods give very similar results. At typical resolutions 
the differences between the two methods at the same resolution are much smaller than between resolutions. The only 
exception from this behavior is at early times, when the iterated method clearly shows first-order convergence that 
goes over smoothly into the expected second-order convergence. 

In the flat empty background, pulses with support well inside x = 1 quickly leave the computational domain. On 
the Choptuik background, we expect to find (exponentially damped) quasiperiodic behavior at late times. We must 
therefore evolve to large values of r (of the order of 10 to 20 background periods A). Not surprisingly we find that 
second-order convergence breaks down after a period or so, both because the quasiperiodic oscillations at different 
resolutions drift out of phase, and because the exponential decay rates are slightly different. At early times, we still 
observe second-order convergence, as described above. 
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